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I. INTRODUCTION 

Monte Carlo simulations, and in particular Markov 
chain based methods, have matured over the last decades 
into a highly versatile and poviferful toolbox for studies 
of systems in statistical and condensed-matter physics 
[U, ranging from classical spin models Q over soft- 
matter problems (3] to quantum systems 0]. Their com- 
petitiveness with other approaches such as, e.g., field- 
theoretic expansions for the study of critical phenomena 
@i 01 1 is IS'igely based on the development and refinement 
of a number of advanced simulation techniques such as 
cluster algorithms Q and generalized-ensemble methods 

Equally important to the generation of simulation 
data, however, is their correct and optimal analysis. In 
this field, a number of important advances over the tech- 
niques used in the early days have been achieved as well. 
These include, e.g., the finite-size scaling (FSS) approach 
PH , turning the limitation of simulational methods to fi- 
nite system sizes into a systematic tool for accessing the 
thermodynamic limit, reweighting techniques [l2| . lifting 
the limitation of numerical techniques to the study of 
single points in parameter space to allow for continuous 
functions of estimates to be studied, as well as advanced 
statistical tools such as the jackknife and other resam- 
pling schemes of data analysis p^ . 

Of these techniques, the statistical data analysis ap- 
pears to have received the least attention. Hence, while 
FSS analyses, even including correction terms, are quite 
standard in computer simulation studies a proper 
analysis and reduction of statistical errors and bias ap- 
pears to be much less common. Here, resampling meth- 



* Electronic address: 
t Electronic address: 



weigel? 



mm-mainz 



janke@itp.um-leipzig.de| 



ods turn out to be very valuable. Although such tech- 
niques offer a number of benefits over more traditional 
approaches of error estimation, their adoption by practi- 
tioners in the field of computer simulations has not yet 
been as universal as desirable. It is our understanding 
that this is, in part, due to a certain lack in broadly ac- 
cessible presentations of the basic ideas which are, in fact, 
very simple and easy to implement in computer codes, as 
is demonstrated below. 

More specifically, data generated by a Monte Carlo 
(MC) simulation arc subject to two types of correlation 
phenomena, namely (a) autocorrelations or temporal cor- 
relations for the case of Markov chain MC (MCMC) sim- 
ulations, which are directly related to the Markovian na- 
ture of the underlying stochastic process and lead to an 
effective reduction of the number of independently sam- 
pled events and (b) cross correlations between different 
estimates extracted from the same set of original time se- 
ries coming about by the origin of estimates in the same 
statistical data pool. The former can be most conve- 
niently taken into account by a determination of the rel- 
evant autocorrelation times and a blocking or binning 
transformation resulting in an effectively uncorrelatcd 
auxiliary time series [l^ . Such analyses are by now stan- 
dard at least in seriously conducted simulational studies. 
On the contrary, the effects of cross correlations have 
been mostly neglected to date (see, however, Refs. [isl - 
ITsj). but are only systematically being discussed follow- 
ing our recent suggestion [l^,[2^. In this article, we show 
how such cross correlations lead to systematically wrong 
estimates of statistical errors of averaged or otherwise 
combined quantities when a nai've analysis is employed, 
and how a statistically correct analysis can be easily 
achieved within the framework of the jackknife method. 
Furthermore, one can even take benefit from the pres- 
ence of such correlation effects for significantly reducing 
the variance of estimates without substantial additional 
effort. We demonstrate the practical relevance of these 
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considerations for a finite-size scaling study of the Ising 
model in two and three dimensions. 

The rest of this article is organized as follows. In Sec. 
II we give a general recipe for a failsafe way of Monte 
Carlo data analysis, taking into account the effects of 
autocorrelations and cross correlations mentioned above. 
After discussing the complications for the more conven- 
tional analysis schemes (but not the jackknife method) 
introduced by histogram reweighting and generalized- 
ensemble simulation techniques in Sec. Ill, we outline the 
role of cross correlations in the process of averaging over 
a set of MC estimates in Sec. IV and discuss the choice 
of an optimal averaging procedure. In Sec. V, these ideas 
are applied to a simulational study of the critical points 
of the two- and three-dimensional Ising models. Finally, 
Sec. VI contains our conclusions. 



In contrast to (O), the estimator O is a random num- 
ber, which only coincides with (O) in the limit TV — >■ oo. 
Under these circumstances, simulational results are only 
meaningful if in addition to the average O we can also 
present an estimate of its variance a^{0). Note that, al- 
though the distribution p{0) of individual measurements 
might be arbitrary, by virtue of the central limit theorem 
the distribution of the averages O must become Gaussian 
for iV — > oo. Hence, the variance cr^(O) is the (only) 
relevant parameter describing the fluctuations of O. If 
subsequent measurements Oi , O2 , • ■ • are uncorrelated, 
we have 
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II. MONTE CARLO ERROR ANALYSIS 

Compared to the task of estimating the uncertainty 
in the result of a lab experiment by simply repeating 
it several times, there are a number of complications in 
correctly determining — and possibly even reducing — 
statistical fluctuations in parameter estimates extracted 
from VICVIC simulations. Firstly, due to the memory 
of the Markovian process, subsequent measurements in 
the time series are correlated, such that the fluctuations 
generically appear smaller than they are. This issue can 
be resolved by a blocking of the original time-series data. 
Secondly, one often needs to know the precision of pa- 
rameter estimates which arc complicated (and sometimes 
non-parametric) functions of the measured observables. 
Such problems are readily solved using resampling tech- 
niques such as the jackknife. 



A. Autocorrelations 

Consider a general Monte Carlo simulation with the 
possible values O of a given observable O appearing ac- 
cording to a probability distribution p{0). This form, of 
course, implies that the system is in thermal equilibrium, 
i.e., that the underlying stochastic process is stationary. 
The probability density p{0) could be identical to the 
Boltzmann distribution of equilibrium thermodynamics 
as for the importance-sampling technique [2l| , but differ- 
ent situations are conceivable as well, see the discussion 
in Sec. IIIII below. If we assume ergodicity of the chain, 
the average 
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for a time series {Oi, O2, . . .} of measurements is an 
unbiased estimator of the mean 



(O) = y dOp(0)0. 
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i.e., ((7^(0)) = CT^(O). This is what we do when estimat- 
ing the statistical fluctuations from a series of indepen- 
dent lab experiments. Markov chain simulations entail 
the presence of temporal correlations, however, such that 
the connected autocorrelation function, 
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is non-zero in general (see, e.g., Ref. 
of the chain implies that Co(s, s + t) = 
Then, the variance of O becomes 
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Monte Carlo correlations decline exponentially, i.e., 

Co(t)'-Co(0)e-*/^"-(O' (5) 

to leading order, defining the exponential autocorrelation 
time Toxp(O). Due to this exponential decay, for N 3> 

t/N of Eq. (gl) from 
and defining the integrated 



Toxp the deviations of the factors 1 
unity can be neglected [! 
autocorrelation time as 
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In view of the 1 /N reduction of variance of the average O 
relative to a single measurement in Eq. ([J) , Eq. ([7]) states 
that the effective number of independent measurements 
in the presence of autocorrelations is reduced by a factor 
of l/2Tint(0). The autocorrelation times Toxp and Tint are 
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FIG. 1: (Color online) Blocking transformation on a time 
series. In the binning analysis, the series is divided into blocks 
of length Nt (a). In the jackknifing analysis, the blocks consist 
of the whole series apart from the entries of a single block (b). 
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not identical, but one can show that the latter is a lower 
bound of the former, Tint(O) < Toxp(O) [25| . 

As long as the autocorrelation time is finite, the dis- 
tribution of averages still becomes Gaussian asymptot- 
ically, such that for N t the variance remains the 
relevant quantity describing fluctuations. To practically 
determine o-^{0) from Eq. an estimate for the auto- 
correlation function is required. This can be found from 
the deflnition ^ by replacing expectation values with 
time averages. It turns out, however, than upon sum- 
ming over the contributions of the autocorrelation func- 
tion for different time lags t in Eq. (j?]) divergent fluctu- 
ations are incurred, enforcing the introduction of a cut- 
off time (26l . [27j . Several approximation schemes have 
been developed using such estimators, but they turn out 
to have severe drawbacks in being computationally ex- 
pensive, hard to automatize and in that estimating their 
statistical accuracy is tedious (see Ref. Q)- 

A more efficient and very intuitive technique for deal- 
ing with autocorrelations results from a blocking trans- 
formation in the spirit of the renormalization group |l4l| 
(in fact, this idea was already formulated by Wilson |28|). 
Much like block spins are defined there, one combines 
Nf, = N/n adjacent entries of the time series. 
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and defines block averages 
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cf. Fig.[IJa). This procedure results in a shorter effective 
time series {O^'' , . . .} with n entries. (We assume 
for simplicity that N is an integer multiple of n.) Obvi- 
ously, the average O and its variance <t^(0) are invariant 
under this transformation. Under the exponential decay 
dS]) of autocorrelations of the original series it is clear (and 
can be shown explicitly however, that subsequent 

block averages O^+i are less correlated than the 

original measurements Ot and Ot+i. Furthermore, the 
remaining correlations must shrink as the block length 
Ni, is increased, such that asymptotically for A''^ — )• oo 
(while still ensuring n ^ 1) an uncorrelated time series 
is produced. Consequently, the nai've estimator ^ can 



FIG. 2: (Color online) Schematic representation of the esti- 
mate (O) of the variance of the average according to Eq. ^ 
for a re-blocked time series as a function of the block length 
Nt. 



be legally used in this limit to determine the variance 
CT^(O) of the average. For the finite time series encoun- 
tered in practice, a block length Nj, ^ t and Ni, N 
must be used. This is illustrated in Figure [2] showing 
the estimate ([2]) for a blocked time series with autocor- 
relation time Tint ~ 13 as a function of the block length 
Nt- It approaches the true variance (J^{0) from below, 
eventually reaching a plateau value where any remaining 
pre-asymptotic deviations become negligible compared to 
statistical fluctuations. If the available time series is long 
enough (as compared to r), it is often sufficient to simply 
lump the data into as few as some hundred blocks and 
restrict the subsequent data analysis to those blocks. As 
a rule of thumb, in practical applications it turns out 
that a time series of length A^ > 10 000 r is required for 
a reliable determination of statistical errors as well as 
autocorrelation times. From Eqs. ^ and (jH]) it follows 
that the integrated autocorrelation time can be estimated 
from 
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within this scheme, where Nt needs to be chosen in the 
plateau regime of Fig. 



B. Covariance and bias 

Apart from providing an estimate of cr^(O) for simple 
quantities, the blocking procedure has the advantage of 
resulting in an effectively uncorrelated auxiliary time se- 
ries which can then be fed into further statistical machin- 
ery, much of which is restricted to the case of independent 
variables. Resampling schemes such as the jackknife [l3| 
provide error and bias estimates also for non-linear func- 
tions of observables without entailing truncation error or 
requiring assumptions about the underlying probability 
distributions. 
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While cr^(O) can be directly computed from the 
blocked time series of O via the estimator ([2]), this ap- 
proach fails for non-linear functions f{{A), (B), . . .) of ex- 
pectation values (A), (B), ... such as, e.g., susceptibil- 
ities or cumulants. A standard approach for such cases 
is the use of error propagation formulas based on Taylor 
expansions p^ . 



estimators (IT3l) as 
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Apart from the truncation error resulting from the re- 
striction to first order in the expansion, this entails a 
number of further problems: if the averages A, B etc. 
are correlated due to their origin in the same simula- 
tion, cross-correlation terms need to be included as well. 
Even worse, for the case of non-parametric parameter 
estimates, such as determining the maximum of some 
quantity by reweighting (see below) or extracting a criti- 
cal exponent with a fitting procedure, error propagation 
cannot be easily used at all. 

Such problems are avoided by methods based on 
repeated sampling from the original data pool, us- 
ing the properties of these meta samples to estimate 
(co-)variance, reduce bias etc. These are modern tech- 
niques of mathematical statistics whose application only 
became feasible with the general availability of comput- 
ers [l^. Most straightforwardly applicable is the jack- 
knife procedure, where the meta samples consist of all 
of the original time series apart from one data block, cf. 
Fig. [TJb). Assume that a set of simulations resulted in a 
collection of time series {O^^i, Ofc,2, • ■ • Ok,Nk}i k = 1, 2, 
. . . for different observables, system sizes, temperatures 
etc. Applying the blocking procedure described above, 
it is straightforward to divide the series in effectively 
uncorrelated blocks. It is often convenient to use the 
same number of blocks n for all series (e.g., 100) which 
can easily be arranged for by the blocking transforma- 
tion as long as Nk / Tk is larger than some minimum value 
(e.g., 10 000) for each simulation and observable. If then 
*Bt = («Bi,t, . . . , *Bfc,t)^ denotes the t^^ block over aU se- 
ries according to Eq. ([8|), where for a constant number 
of blocks the block lengths iV^ fc = Nk/n might vary be- 
tween the different scries under consideration, one defines 
the corresponding jackknife block as the complement 
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cf. Fig. [TJ Considering now an estimator 9{{Ot}) for 
some parameter 9 depending on (some or all of) the dif- 
ferent series, we define the corresponding estimates re- 
stricted to jackknife block -3^. 



(13) 



The variation between these estimates taken from the 
same original data can be used to infer the sample vari- 
ance. If one denotes the average of the jackknife block 
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an estimate for the sample variance of the estimator Q is 
given by [29| 
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This is very similar to the simple estimate ^ for the vari- 
ance of the average, but it comes with a different pref- 
actor which serves a twofold purpose: it reweights the 
result from the effective jackknife series of length n — 1 
to the original length n and takes care of the fact that all 
of the jackknife block estimates 0(s) are strongly corre- 
lated due to them being based on (almost) the same data. 
The general Eq. (|15p forms a conservative and at most 
weakly biased estimate of the true variance [l^, which 
lacks the truncation error of schemes based on Eq. (fTTj) 
and is applicable to non-parametric parameter estimates. 

In a slight generalization of Eq. ([15]) it is possible to 
also estimate covariances. For a number of estimators 
di{{Ot}), i = 1, 2, . . ., a robust jackknife estimator of 
the covariance matrix is given by 
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In a similar way the bias of estimators can be reduced, 
i.e., deviations between the mean of an observable and 
the expectation value of some estimator that disappear 
with increasing sample length. For a detailed discussion 
we refer the reader to Refs. [isl. [29j. 

A general procedure for the analysis of simulation data 
based on blocking and jackknife techniques hence has the 
following form: 

1. Decide on the number n of jackknife blocks to be 
used. For most purposes, of the order of 100 — 500 
blocks are sufficient. 

2. For each original time series recorded in a collec- 
tion of simulations, examine the block averages © 
as a function of the block length N/n. If the re- 
sult for n blocks is in the plateau regime of Fig. [5] 
everything is fine; otherwise, one needs to record a 
longer time series (and possibly take measurements 
less frequently to keep the amount of data manage- 
able) . 

3. For each parameter to be estimated, compute the n 
jackknife block estimates ((T5)) as well as the average 
(|14l) and combine them to calculate the variance 
(|15p . For a number of different parameter estimates 
9i, the jackknife block estimates can also be used 
to calculate the covariance (HI 
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III. HISTOGRAMS AND ERRORS 

An increasing number of successful Monte Carlo tech- 
niques rely on reweighting and the use of histograms 
[l|. This includes the (multi-) histogram method of 
Refs. [HillOl as well as the plethora of generalized ensem- 
ble techniques ranging from multicanonical simulations 
to Wang-Landau sampling [l^l- Such methods arc 
based on the fact that samples taken from a known prob- 
ability distribution can always be translated into samples 
from another distribution over the same state space. As- 
sume, for simplicity, that states are labeled {si} as ap- 
propriate for a spin system. If a sequence {si}t, t = 1, 2, 
. . . was sampled from a stationary simulation with proba- 
bility density Psim({si}), an estimator for the expectation 
value of the observable O relative to the equilibrium dis- 
tribution is given by 



N 



o = 



t=l 



Pcq{{Si}t) 
Psimi{si}t) 



N 

E 

t=i 



Peq{{Si}t) 
Psim{{Si}t) 



(17) 



For a finite simulation this works as long as the sampled 
and the equilibrium distributions have sufficient overlap, 
such that the sampled configurations can be represen- 
tative of the equilibrium average at hand. For simple 
sampling one has Psim = const and hence must weight 
the resulting time series with the Boltzmann factor 



Poq{{Sl}) 



(18) 



where '^({si}) denotes the energy of the configuration 
{si} and is the partition function at inverse tempera- 
ture /3 = l/kBT. For importance sampling, on the other 
hand, Psim = Peq, such that averages of time series are di- 
rect estimates of thermal expectation values. If samples 
from an importance sampling simulation with psim = PPo 
should be used to estimate parameters of Peq = Pp, 
Eq. (|17p yields the familiar (temperature) reweighting 
relation 



Ofi = 



(19) 



where Et = H({si}t). Completely analogous equations 
can be written down, of course, for reweighting in param- 
eters other than temperature. Similarly, canonical aver- 
ages at inverse temperature /3 arc recovered from mul- 
ticanonical simulations via using Eq. (jl7p with Pshn = 

Pmuca and Peq = P/S- 

Reliable error estimation (as well as bias reduction, 
covariance estimates etc.) for reweighted quantities is 
rather tedious with traditional statistical techniques such 
as error propagation (3l1 |. Resampling methods, on the 
other hand, allow for a very straightforward and reliable 
way of tackling such problems |52| . For the jackknife 



approach, for instance, one computes jackknife block es- 
timates of the type (|17p by simply restricting the set of 
time series to the s"^ jackknife block -3 s. With the jack- 
knife average ((T4| . e.g., the variance estimate ((T5|) with 
9 = O can be straightforwardly computed. Similar con- 
siderations app ly to covariance estimates or bias reduced 
estimators [13| . Extremal values of thermal averages can 
be determined to high precision from the continuous fam- 
ily of estimates where error estimates again follow 
straightforwardly from the jackknife prescription. 



IV. VARIANCE REDUCTION 

Temporal correlations resulting from the Markovian 
nature of the sampling process have been discussed in 
Sec. Ill Al above, and we assume that they have been ef- 
fectively eliminated by an appropriate binning procedure. 
Extracting a number of different parameter estimates 9i, 
1 = 1,2,... from the same number of original simulations 
it is clear, however, that also significant cross correlations 
between estimates 9i and 9j can occur. These have pro- 
found consequences for estimating statistical error and 
reducing it by making the best use of the available data 

If a given parameter estimate 9 depends on several ob- 
servablcs of the underlying time series that exhibit cross 
correlations, this fact is automatically taken into account 
correctly by the jackknife error estimate p^. This is in 
contrast to error analysis schemes based on error propa- 
gation formulae of the type ([TT|). where any cross corre- 
lations must be taken into account explicitly. Insofar the 
outlined approach of data analysis is failsafe. We want to 
go beyond that, however, in trying to optimize statistical 
precision of estimates from the available data. If we at- 
tempt to estimate a parameter 9, we ought to construct 
an estimator 

e^T{{dt}), 

which is a function of the underlying time series with the 
property that {9) = 9 (at least for n — > oo). Obviously, 
there usually will be a large number of such functions T 
and it is not possible, in general, to find the estimator 
9 of minimal variance. We therefore concentrate on the 
tractable case where is a linear combination of other 
estimators 9i, i = 1, . . . , 
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There are different possibilities to ensure the condition 

{9) = 9: 

1. All estimators have the same expectation, {9i) = 9, 
and = 1- 
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2. One estimator is singled out, say (6*1) = 9, ai = 1, 
and the rest has vanishing expectation, (di) ~ 0, ai 
arbitrary, i > 2. 

3. More comphcatcd situations. 

The first type describes the case that we have several 
different estimators for the same quantity and want to 
take an average of minimum variance [l9[. The second 
case is tailored for situations where existing symmetries 
allow to construct estimators with vanishing expectation 
whose cross correlations might reduce variance [20l |. 

To optimize the analysis, the parameters ai in (j20p 
should be chosen such as to minimize the variance 

fc k 
ij = l i.j = l 

(21) 

For case one above, we introduce a Lagrange multiplier to 
enforce the constraint ai = 1, and the optimal choice 
of ai is readily obtained as 



leading to a minimum variance of 

'^'W = • 

Very similarly, case two leads to the choice [SS 

k 

a, = -Y,[V'{e)-%,m,,, 

j=2 



(22) 



(23) 



(24) 



where r'(^) denotes the submatrix of \^{0)]ij with i, j > 
2. Since the formalism for both cases is practically iden- 
tical, in the following we will concentrate on case one. 

Let us take the time to compare the optimal choice of 
weights expressed in Eqs. ([^ and ([21]) with that used 
in more traditional approaches. Ignoring the presence of 
cross correlations, several parameter estimates are often 
combined using an error weighting scheme only, i.e., by 
choosing weights 



(25) 



While the more general expression ((22|) reduces to the 
weights (PS)) in the absence of correlations, the choice (PS)) 
is not optimal as soon as cross correlations are present. 
Still, the resulting average 9 remains a valid estimator of 
the parameter 0. In contrast, the usually used variance 
estimate derived from the expression 
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is no longer even correct when cross correlations come 
into play. As will be seen below from the discussion of 
Ising model simulations in Sec. |Vl cfI^^^o-cA^) generically 
leads to underestimation of the true variance, but occa- 
sionally over-estimates are possible as well. 

The practical implementation of the described scheme 
of weighting and error analysis is straightforward with 
the toolset outlined in the previous sections. The co- 
variance matrix of the estimates 9i is readily computed 
via the jackknife expression (|16|) . This allows to estimate 
the optimal weights from inserting in Eq. (|22p (or the 
analogue ([M)) for case two) and the variance of the result- 
ing optimal estimator is determined from the expression 
(|23p . In total, the necessary analysis can be summarized 
as follows: 

1. Perform a binning analysis to see whether for a 
given number of blocks n the block averages © for 
all time series at hand are effectively uncorrelated. 

2. For each parameter estimate 9i compute the n jack- 
knife block estimates as well as their average 
and estimate their covariance matrix from Eq. (jl6|) . 

3. For those estimates Oi to be combined into an av- 
erage ^, an estimate of the optimal weighting pa- 
rameters ai is given by Eq. (|22p with the estimate 
Vij calculated in the previous step. Likewise, the 
variance of the resulting average is estimated from 
Eq. (IMD. 

In some cases, it is necessary to already have variance 
estimates of intermediate data available for properly de- 
termining the jackknife block estimates This typ- 
ically occurs when 9i is a parameter resulting from an 
(ideally error weighted) fit to a number of data points, 
such as for the case of a critical exponent, see the discus- 
sion below in Sec. |Vl In these cases it is straightforward 
to iterate the jackknifing procedure to second order by 
considering each jackknife block as the initial time series 
of another jackknife analysis (32| . 

In view of the sometimes counter-intuitive results of 
computing weighted averages taking cross correlations 
into account (see the results in Sec. fVl below), it is in- 
structive to examine the simple case of just two different 
estimates 9i and 6*2 in somewhat more detail. This is 
done in Appendix VK\ 



V. APPLICATION TO THE ISING MODEL 

Although the outlined scheme of Monte Carlo data 
analysis is completely general, it is useful to see how 
it works out for a specific example. In particular, one 
would like to know if the typical cross correlations are 
sufficiently strong to have significant impact on the re- 
sults. To answer this question, we performed a finite-size 
scaling (ESS) analysis of the ordering transition of the 
ferromagnetic Ising model in two and three dimensions. 
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FIG. 3: (Color online) Fits of the functional forms pO|) 
resp. (|3ip to the maxima data ^max of the following quan- 
tities computed for the case of the 2D Ising model: A — 

dln^^ dln^^ dl^^ ^ ^ (f^^^ 

tom). For performing the fits, the correction terms in brack- 
ets of Eqs. (|30p and (|31[l were neglected. The actual fits have 
been performed on the size range 32 < L < 192. The slopes 
of the lines are identical to the inverse of the corresponding 
estimates of the correlation length exponent f listed in Ta- 
ble H 



A. Simulation details 

We studied the critical behavior of the nearest- 
neighbor, zero-field, ferromagnetic Ising model with 
Hamiltonian 



= ±1 



(27) 



on square and simple cubic lattices of edge length L, us- 
ing periodic boundary conditions. Close to criticality, an 
importance-sampling Monte Carlo simulation with local 
update rule suffers from critical slowing down, r ~ L^, 
with a dynamical critical exponent z « 2. To alleviate 
this problem, we used the single-cluster update algorithm 
[33| resulting in a dramatic speed-up of the relaxation 
process. For two and three dimensions we performed sim- 
ulations at a fixed temperature close to the asymptotic 
critical temperature for a number of different system sizes 
to enable a systematic FSS study. The raw data consisted 
of time series with 4 x 10^ approximately independent 
samples of the configurational energy and magnetization 
for each system size under consideration. Using the jack- 
knifing analysis described above, these original time se- 
ries were then analyzed using n effectively uncorrelated 
bins, where n — 100 was chosen unless stated otherwise. 



B. Finite-size scaling analysis 

There is now a broad consensus that finite-size effects 
are (in most cases) not merely a drawback of approaches 
depending on finite system sizes, but can be turned into a 



powerful tool for extracting the asymptotic behavior [11 
A number of different practical implementations of this 
idea in terms of specific FSS schemes have been derived 
and successfully applied to the analysis of critical phe- 
nomena, see, e.g., Refs. [33 - [36j . Although our consider- 
ations regarding the data analysis apply rather generally 
to all these techniques, for illustrative purposes we con- 
centrate here on the rather popular method outlined in 
Ref. [13] . It is focused on the analysis of the locations and 
values of extrema of standard thermodynamic quantities 
such as the specific heat, magnetic susceptibility, cumu- 
lants etc. According to the theory of finite-size scaling 
[O, [13 ' locations of such pseudo-critical points are 
shifted away from the true critical coupling /3c according 
to 

/3(^max, L)=I3, + AnL-\l + A,L-^ + •••), (28) 

where A denotes an observable with a pseudo-critical 
maximum such as the specific heat (for a > 0). The 
generic value for the shift exponent A predicted by FSS 
theory is A = where v is the correlation length expo- 
nent (for exceptions see, e.g., Ref. [13 )• A drawback of 
using Eq. ([^5]) directly is that non-linear fits in the three 
parameters /3c, a-nd A resp. v are required even for the 
simplest case of ignoring the correction-to-scaling terms 
in brackets. To alleviate this problem, it has been sug- 
gested to consider quantities such as the magnetization 
cumulants [Hi 



1 



(|m 



2i\ 



3(|r 



1,2,3, 



(29) 



for which the maxima of the temperature derivatives have 
a critical scaling form 



dU2^ 



d/3 



C/,.oi^/''(l + C/.,ci-" + •••), (30) 



and hence allow to determine v without prior knowledge 
of the transition coupling /3c. If, again, the correction 
terms in brackets are ignored, this form even represents 
a linear fit (in logarithmic representation) resulting in 
very stable results. While initially only the fourth-order 
cumulant ?74 was considered, the authors of Ref. [13| sug- 
gested to use a variety of different cumulants U2i with 
I = 1, 2, 3, ... to improve the accuracy of the v es- 
timate. A number of further quantities with the same 
scaling behavior can be constructed, for instance loga- 
rithmic temperature derivatives of the magnetization. 



dln(|TO|' 



d/3 



A,oi'/''(i + A,ci" 



(31) 



which for i = 1, 2, . . . yields another series of v estimates. 

Once v has been determined from fits of the functional 
forms (|30p and ([3T|) to the data, one might return to the 
shift relation ([25)1 and (assuming A = 1/;^) determine the 
transition coupling /3c from linear fits with a fixed value 
of V. Finally, the remaining standard critical exponents 
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TABLE I: Fit parameters and correlation data for estimating the critical exponent v from single-cluster update Monte Carlo 
simulations of the 2D Ising model. The exponent estimates are extracted from fits of the functional forms (|30p and (|31|l to the 
data, neglecting the correction terms in the brackets. Deviations from the exact value v = \ are computed relative to = 1 
(Arci) and in multiples of the estimated errors listed in the column labeled "a" (A^). 



fits correlation coefScients/weights 





^min 




V 


a 


Arcl 


A,, 


Q 


d.o.f. 


dln( m ) 
d/3 


dln(m^) 
d/3 


dln( mp) 

d/3 


^ 


dC/4 
^ 


H ln/iml\ 

dp 
dln(|m|3) 


32 


192 


1.0085 


0.0183 


0.85% 


0.47 


0.52 


4 


1.0000 


0.9743 


0.9385 


0.9197 


0.8971 


32 
32 


192 
192 


1.0128 
1.0175 


0.0194 
0.0201 


1.28% 
1.75% 


0.66 
0.87 


0.47 
0.40 


4 
4 


0.9743 
0.9385 


1.0000 
0.9910 


0.9910 
1.0000 


0.8167 
0.7431 


0.8687 
0.8198 


difa 
d/3 


32 


192 


1.0098 


0.0281 


0.98% 


0.35 


0.57 


4 


0.9197 


0.8167 


0.7431 


1.0000 


0.8596 


dC/4 
d/3 


32 


192 


1.0149 


0.0511 


1.49% 


0.29 


0.70 


4 


0.8971 


0.8687 


0.8198 


0.8596 


1.0000 


^plain 




f^uncorr 
^corr 


1.0127 


0.0141 
0.0269 


1.27% 
1.27% 


0.90 
0.47 






1.0000 


1.0000 


1.0000 


1.0000 


1.0000 






'^uncorr 
^^corr 


1.0123 


0.0102 
0.0208 


1.23% 
1.23% 


1.21 
0.59 






0.3145 


0.2714 


0.2483 


0.1322" 


0.0336 






fJcorr 


0.9935 


0.0078 


-0.65% 


-0.84 






5.0067 


-2.4259 


-0.2807 


-1.1958 


-0.1043 



"Note that the similar Tabic II of Ref. [l9l | contains a mistake in the last two lines, where the weights 0.1322 and 0.0336 as well as 
— 1.1958 and —0.1043 appear interchanged with respect to the correct data represented here. 



can be estimated from the well-known FSS forms of the 
specific heat cy, the magnetization m and the magnetic 
susceptibility x, 



Cy|„.ax=CoL"/^(l+C,i-"' + ...), 

(|m|)inf = 7noi-^/''(l + muL-^ + •••), 

Xl„,ax = X0i^/^(l+Xfci'"' + •••), 



(32) 



where (|m|)inf denotes the (modulus of the) magnetiza- 
tion at its inflection point. The directly estimated expo- 
nents arc therefore v and the FSS exponents ajv^ P/v 
and ^/v, which can be combined to yield a, /3, 7 and 
v. The remaining exponents 5 and 77 are not directly de- 
termined here; instead we assume that their values are 
deduced from the exponents a, /?, 7 and v via standard 
scaling relations. 

For determining the (location and value of the) max- 
ima occurring in Eqs. (gHl, O, (EH) and we used 
the rcweighting technique outlined in Sec. IIIII starting 
from the data of a single simulation per system size per- 
formed at or close to the asymptotic transition point. 
The derivatives with respect to /3 in Eqs. (pO)) and ((3T|) 
are easily shown to be equivalent to combinations of mo- 
ments of energy and magnetization at a single, fixed tem- 
perature (3^ . For the fourth-order cumulant L/4, for in- 
stance, one has 



dC4 
"d^ 



2(m'*)[(m-^)(e) — {m?e)] — {m?)[{m'^) {e) — (rn^e)] 



3(m2) 



2\3 



(33) 

Therefore no numerical differentiation is required when 
making use of the relations (PU)) and (PT|) . In some sit- 
uations it might be impractical to store the whole time 



series of original measurements of internal energy and 
magnetization, in which case the exact reweighting rela- 
tion (|19p might be replaced by a Taylor expansion with 
respect to j3 around the simulation coupling /3o, where 
then cumulants of e and m appear as the expansion co- 
efficients. In most cases, however, it is much simpler and 
more versatile in terms of the data analysis to work with 
the original time series. In view of the typically available 
storage resources today, this approach should be com- 
fortably feasible in most situations. 



Considering the set of critical exponents a, /3, 7, v 
(as well as ry and <5), it is useful to recall that they are 
subject to a number of exact scaling relations, namely 
the Rushbrooke identity a + 2/3 -I- 7 = 2, Fisher's scaling 
law 7 = v{2 — rf), the relation a + /3(1 + 5) — 2, as well 
as (in most cases) the hyperscaliiig relation a — 2 — Av^ 
where d is the spatial dimension [40|. As a consequence 
of these four equations, only two of the six original ex- 
ponents are independent. While the scaling relations 
have been occasionally used to check a set of indepen- 
dently estimated exponents for consistency, we would like 
to point out that the existence of these exact relations 
should rather be used for improving the ■precision of ex- 
ponent estimates. In particular, in the language of the 
renormalization group, it is natural to express the con- 
ventional scaling exponents in terms of the scaling dimen- 
sions Xt and Xfi of the operators coupling to temperature 
and magnetic field |4l| . respectively, which are the only 
relevant operators for the Ising model [43 |. For the four 
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TABLE II: Fitting and averaging results for estimating the critical coupling /3c of the 2D Ising model from the shifts of pseudo- 
critical temperatures according to Eq. (|28p . For performing the fits, the correlation length exponent was fixed at its exact value 
v = 1. The column Aroi indicates the relative deviation of the estimates from the exact result Pc — ^ ln(l + y/2) « 0.4406868. 





Po 


fits 
a 




Q 


Cv 


d(l,n|) 
dp 


correlation coefficients/ weij 
dln(jm|) dln{m2) dln{m^) 
dp dp dp 


j;hts 
df/2 
dp 


dC/4 

"d^ 


X 


Cv 


0.440709 


0.000101 


0.0051% 


0.35 


1.0000 


0.7881 


0.3965 


0.3645 


0.3569 


0.3502 


0.2599 


0.1394 


d(|mj) 
d/? 


0.440798 


0.000073 


0.0251% 


0.09 


0.7881 


1.0000 


0.7116 


0.6417 


0.6043 


0.7166 


0.5345 


0.6236 


d ln(|m|) 

dp 
dln(m2) 

dp 
dln(m'^) 


0.440711 


0.000408 


0.0055% 


0.56 


0.3965 


0.7116 


1.0000 


0.9740 


0.9365 


0.9122 


0.8613 


0.6477 


0.440799 


0.000504 


0.0254% 


0.42 


0.3645 


0.6417 


0.9740 


1.0000 


0.9899 


0.8119 


0.8111 


0.5264 


0.440918 


0.000567 


0.0525% 


0.29 


0.3569 


0.6043 


0.9365 


0.9899 


1.0000 


0.7415 


0.7608 


0.4554 


dSi 


0.440576 


0.000338 


-0.0251% 


0.70 


0.3502 


0.7166 


0.9122 


0.8119 


0.7415 


1.0000 


0.8854 


0.8413 


dp 

X 


0.440308 


0.000708 


-0.0860% 


0.94 


0.2599 


0.5345 


0.8613 


0.8111 


0.7608 


0.8854 


1.0000 


0.6265 


0.440699 


0.000045 


0.0028% 


0.67 


0.1394 


0.6236 


0.6477 


0.5264 


0.4554 


0.8413 


0.6265 


1.0000 



/Sc.piain (Tuncorr 0.440690 0.000151 0.0007% 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 

(Tcorr 0.000322 0.0007% 

Pc.orr (Tuncorr 0.440725 0.000036 0.0086% 0.1364 0.2530 0.0073 0.0049 0.0039 0.0105 0.0024 0.5815 

CTcorr 0.000059 0.0086% 

/3c,cov o-corr 0.440687 0.000018 0.0001% 0.2601 -0.3347 0.2086 -0.0566 -0.0186 -0.2997 0.0276 1.2133 



exponents considered here, this means that 

— = d - 2xt, - = Xh, 

\ (34) 

- = d - 2xh, - = d- xt, 
u u 

such that Xh can be independently estimated from x^ = 
15 /v and Xh = d/2 — 7/2j/, whereas xt might be deter- 
mined from Xt ~ d — 1 / V as well as Xt = d/2 — a/2v. 



C. Two-dimensional Ising model 

For the two-dimensional (2D) Ising model, single- 
cluster update simulations were performed for square lat- 
tices of size L = 16, 24, 32, 48, 64, 96, 128 and 192. All 
simulations were done directly at the asymptotic critical 
coupling /3c = 5 ln(l + ^2) « 0.4406868. For the range 
of system sizes under consideration, it turned out that 
simulations at this single temperature were sufficient for 
reliably studying the pseudo-critical points defined by the 
maxima of the various quantities under consideration by 
reweighting, i.e., the overlap of histograms between the 
simulation and analysis temperatures turned out to be 
sufficiently large. 

We first extracted a number of estimates of the corre- 
lation length exponent v from investigating the maxima 
of logarithmic derivatives of magnetization moments for 
i = 1, 2 and 3 as well as the maxima of the derivatives 
of the second-order and fourth-order cumulants U2 and 
U4 using the reweighting scheme outlined above. The 



locations and values of the maxima themselves were de- 
termined by a golden section search algorithm j43| . The 
resulting maxima as a function of system size are shown 
in Fig. El together with fits of the forms and ^ 
to the data. Here, we used fits without the correction 
terms in the brackets of Eqs. ([501 and (|5T|) on the fit 
range L > 32, which works very well as is apparent from 
the presentation in Fig. [3] and the corresponding values of 
the quality-of-fit parameter Q listed in the eighth col- 
umn of Table m The fourth and fifth column contain the 
resulting estimates of the exponent u together with the 
statistical errors estimated from a weighted least-squares 
fitting procedure [i^ . A glance at Table H] reveals that 
all single estimates are statistically consistent with the 
exact result 1/ = 1, but they exhibit a rather large vari- 
ation in statistical accuracy with the biggest statistical 
error being almost three times larger than the smallest. 
We use the jackknife estimator (fT6|) and a second-order 
jackknifing procedure to estimate the statistical correla- 
tions of the individual estimates of ly. The data on the 
right hand side of Table [J showing the correlation coef- 
ficients p ~ Tij/aiUj for the different estimates reveal 
that correlations between all pairs of estimates are large 
with p > 0.8. With all estimates being derived from simi- 
lar expressions containing magnetic moments, this result 
probably does not come as a surprise. 

Under these circumstances, one might wonder whether 
it is worthwhile to attempt a linear combination of the 
form (|20[) of the various v estimates rather than quoting 
the single most precise estimate as final result, which in 
the present case is given by the value v = 1.0085(183) 
resulting from the FSS of dln(|TO|)/d/3. For the purpose 
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TABLE III: Averaging results for estimates of the critical 
exponent v and the critical coupling /3c of the 2D Ising model 
from non-linear three-parameter fits to of the functional form 
(|28|) . The observables used are those listed in Table [III 







V 


a 


/3c 


a 


^plain 


f^uncorr 


0.8101 


0.0428 


0.439491 


0.000337 




C^corr 




0.0973 




0.000715 


^crr 


^uncorr 


0.8949 


0.0228 


0.440295 


0.000099 




f^corr 




0.0435 




0.000169 


^cov 


(^corr 


0.9980 


0.0148 


0.440658 


0.000072 


exact 




1.0000 




0.440687 





of combining estimates, we consider the traditional ap- 
proaches of taking a "plain average i^piain with 

«f = \ (35) 

as well as the error-weighted average of Eq. (|25p and 
compare them to the truly optimal covariance-weighted 
average t'cov defined by the weights of Eq. ([22]) . Ignoring 
the presence of correlations (as was the case in most pre- 
vious studies), one would estimate the error associated 
to the plain average as 

i 

and, likewise, the variance of the error- weighted average 
is given by (Junmrr defined in Eq. p6|. The true vari- 
ances of t'piain and Uerr in the presence of correlations, on 
the other hand, can also be easily derived formally, and 
will contain the elements of the covariance matrix F of 
the individual estimates I'i. From the practical perspec- 
tive, the jackknifing analysis outlined here automatically 
takes those correlations into account. We refer to these 
correctly defined variances with the notation ct^q^j.. 

The plain, error-weighted and covariance-weighted av- 
erages for 1/ with the corresponding variance estimates 
are listed in the lower part of columns four and five of 
Table [D As with the individual estimates, each of the 
three averages is statistically compatible with the exact 
result V = 1. While the nai've error estimates CTuncorr seem 
to indicate that performing the plain or error-weighted 
average reduces statistical fluctuations compared to the 
single estimates, taking correlations into account with 
the jackknifing scheme resulting in cr^^^^ reveals that vari- 
ances are grossly underestimated by (J^^^on and, in fact, 
compared to both the plain (ctcoit = 0.0269) and error- 
weighted ((Tcorr = 0.0208) averages the single estimate 
of V stemming from dln(|m|)/d/3 has smaller statistical 
fluctuations (cr = 0.0183). Performing those averages 
therefore decreases precision instead of improving it! The 
truly optimal average of Eq. ((22|) . on the other hand, 
results in the estimate i/ = 0.9935(78), whose fluctua- 
tion is about 2-3 times smaller than those of the error- 
weighted average and the single most precise estimate. 



The reduced variance of this last estimate seems to be 
corroborated by the smallest deviation also from the ex- 
act result I' = 1. A glance at the data collected in Table 
U reveals that, somewhat astonishingly, the optimal av- 
erage is smaller than all of the individual estimates of i/ 
(see also the graphical representation of this fact in Fig. 1 
of Ref. [l^). This situation which, of course, can never 
occur for the error-weighted average where all weights 
< ctj < 1, is connected to the fact that the more gen- 
eral weights of Eq. ()22[) are unbounded and, in particular, 
can become negative. This fact reflects in the computed 
weights for the different averaging schemes collected in 
the lower right hand part of Table ID Clearly, the weights 
for the error-weighted and covariance-weighted averages 
are dramatically different and, in particular, some of the 
latter turn out to be negative. It is intuitively clear that 
such negative weights are necessary to cancel the effects 
of strong mutual correlations. The asymmetry in the 
weights leading to the possibility of the average lying 
outside of the range of the individual estimates results 
from the asymmetry of the individual variances in con- 
nection with the cross correlations. This effect can be 
explicitly understood for the case of only two estimates, 
cf. the discussion in Appendix Rl 

One might wonder whether the suggested weighting 
scheme requiring to estimate the full covariance matrix 
is statistically robust. It is clear, for instance, that the 
jackknife estimator (|16|) for the covariance will become 
more precise as more jackknife blocks are used — at 
the expense of an increased computational effort. To 
check for such effects we repeated our analysis while us- 
ing n = 200 instead of n = 100 jackknife blocks. Most of 
the estimates for the correlation coefficients are almost 
unchanged by this new analysis with the largest devia- 
tion being of the order of 3%. The same holds true for 
the resulting weights in the optimal average, where only 
the weight of the estimate resulting from d ln(|7Tip)/d/3 
changes substantially from a = —0.2807 to a = —0.5499. 
The final optimal estimate D ~ 0.9908(78) is fully com- 
patible statistically with the analysis using 100 jackknife 
blocks. Using (as a consistency check) completely in- 
dependent simulations for producing the individual esti- 
mates of ly, on the other hand, indeed results in a unit 
matrix of correlation coefficients within statistical errors 
and, consequently, the error-weighted and covariance- 
weighted averages coincide in this limit. Finally, we also 
find that the numerical inversion of the covariance matrix 
required for computing the weights in Eq. is in gen- 
eral stable and unproblematic. It is clear, however, that 
in the presence of very strong correlations the resulting 
weights of individual estimates will depend sensitively 
on the entries of the covariance matrix, since in the limit 
of perfect correlations all choices of weights become de- 
generate, see also the discussion of the case of only two 
estimates in Appendix lAl 

We now turn to the determination of the transition 
coupling f3c from the shift relation ([28]) . We considered 
the locations of the extrema of the specific heat cy, the 
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slope d(|m|)/d/3 of the (modulus of the) magnetization, 
the logarithmic derivatives dln(|mp)/d/3 for i = 1, 2 and 
3, the cumulant derivatives d?72/d/3 and dUi/dfi as well 
as the magnetic susceptibility x- In order to most clearly 
demonstrate the effects of the present correlations, we 
first performed fits of the form ([SS)) using the exact cor- 
relation length exponent v ~ 1. The corresponding fit 
results are collected in Table HIl For the fits wc ignored 
the correction terms indicated in the brackets of Eq. p8| , 
leaving out the smallest system sizes instead. This ap- 
proach appears justified in view of the good fit qualities 
refiected in the Q values of Table HIl As for the fits for de- 
termining V, all single estimates of /3c are consistent with 
the true asymptotic values of /3c ~ 0.4406868 within error 
bars. The corresponding standard deviations, however, 
vary dramatically, decreasing by a factor of 15 from the 
estimate resulting from dC/4/d/3 to that of the suscepti- 
bility X- The results of the correlation analysis are pre- 
sented on the right-hand side of Table while the log- 
arithmic magnetization derivatives and cumulants again 
show very strong correlations, the results of the remain- 
ing quantities are somewhat more independent, showing, 
in particular, a rather clear separation of the energetic 
from the magnetic sector. For the averages of single esti- 
mates the present correlations again lead to a significant 
underestimation of the true variance for the plain and 
error- weighted cases and, in fact, both of them are less 
precise than the best single estimate stemming from the 
scaling of the susceptibility x, cf. the data in the lower 
part of Table |TT1 The truly optimal average of Eq. ((^^ 
results in /3c = 0.440687(18), where the statistical error is 
about threefold reduced compared to the error-weighting 
scheme. Very similar results are found when using the 
value V = 0.9935(78) found from the analysis summarized 
in Table U where we arrive at a final covariance- weighted 
average of /3c = 0.440658(17) [35]. Here, the second er- 
ror estimate in square brackets refers to the sensitivity of 
the result for /3c to the uncertainty in v indicated above, 
which turns out to be symmetric with respect to upwards 
and downwards deviations of v here. 

As an alternative to the two-step process of first deter- 
mining 1/ from the relations (|30l) and ([31]) and only after- 
wards estimating /3c from Eq. ()28p . one might consider 
direct fits of the form ((25)) to the maxima data of the 8 
observablcs listed above determining v and /3c in one go. 
Here, again, fits on the range 32 < L < 192 neglecting 
any corrections to the leading scaling behavior are found 
to be sufficient. The results for the plain, error- weighted 
and covariance- weighted averages for both parameters, v 
and /3c, are collected in Table IIIII Consistent with the 
previous results, it is seen that neglecting correlations in 
error estimation leads to a sizable underestimation of er- 
rors and, on the other hand, using the optimal weighting 
scheme of Eq. ([22]) statistical errors are significantly re- 
duced, an effect which is also nicely illustrated by the 
very good fit of the resulting parameter estimates with 
the exact values. 

Finally, we turn to the determination of the remain- 



TABLE IV: Determining the magnetic and energetic scal- 
ing dimensions Xh and Xt of the 2D Ising model by weighted 
averages over various individual estimates. 







Xh 


a 


Xt 


a 




f^uncorr 


0.1219 


0.0027 


1.0085 


0.0117 




f^corr 




0.0021 




0.0213 


Ocrr 


f^uncorr 


0.1261 


0.0016 


1.0048 


0.0082 




f^corr 




0.0013 




0.0136 




f^corr 


0.1250 


0.0010 


1.0030 


0.0096 


exact 




0.1250 




1.0000 





ing critical exponents. As outlined above, wc do this 
by combining different estimates using covariance anal- 
ysis to improve the results for the scaling dimensions, 
thus ensuring that the scaling relations arc fulfilled ex- 
actly. From a glance at Eq. ([M)) one reads off that the 
magnetic scaling dimension Xh can be determined from 
Xh ~ P/v and Xh = d/2 — 7/2z^. We therefore determine 
P/v from the ESS of the (modulus of the) magnetization 
at its infiection point and estimate ^/v from the ESS of 
the susceptibility maxima, resulting va j3/v ^ 0.1167(54) 
and ^/v = 1.7458(40), respectively. As the correlation 
analysis reveals, the two resulting estimates of Xh are 
anti-correlated to a considerable degree with correlation 
coefficient —0.64. As a consequence, conventional error 
analysis neglecting correlations ower-estimates statistical 
fiuctuations. Still, choosing optimal weights according to 
Eq. (P^ is able to reduce variance, resulting in a com- 
bined estimate Xh = 0.1250(10) right on top of the exact 
result Xh = 1/8, cf. the data collected in Table HVl The 
energetic scaling dimension xt, on the other hand, might 
be computed from Xt = c?— 1 /z^ as well as Xt = djl—ajlv . 
We therefore use the five individual estimates of v listed 
in Table |T] as well as the FSS of the maximum of the 
specific heat to estimate xt . The latter fits are somewhat 
problematic due to the logarithmic singularity of the spe- 
cific heat corresponding to ajv = 0, and it turns out that 
a fit of the form 

cy.max = Co + ciL"'/'' In L 

including a scaling correction is necessary to describe the 
data. Combining all individual estimates in an optimal 
way, we arrive at Xt = 1.0030(96), well in agreement with 
the exact result Xt = 1, cf. the right hand side of Table 

CYl 



D. Three-dimensional Ising model 

Cluster-update simulations of the ferromagnetic Ising 
model in three dimensions (3D) were performed for 
simple-cubic lattices of edge lengths L = 8, 12, 16, 24, 
32, 48, 64, 96 and 128. All simulations were performed at 
the coupling /3 = 0.221 654 9 reported in a high-precision 
study as estimate for the transition point |45| . since it 
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TABLE V: Fit parameters and correlation data for estimating the critical exponent v from single-cluster update Monte Carlo 
simulations of the 3D Ising model. Fits of the functional form (|3ip including the correction term were used for the logarithmic 
magnetization derivatives d ln{|mP)/d/3 for i = 1, 2 and 3, while fits of the form (|30|l without correction term were used for the 
derivatives of the cumulants f/2 and Ui_. The relevant reference values vs, v = 0.6301(4) taken from Ref. [4^ . 



fits 



d.o.f. 



dln(|m|) 

dln(rn^) 

d/3 
dln{|m|^ 

difa 
d/3 
dCA 

d/3 



correlation coefficients/weights 
dln(|mi) dln(m^) dln(|?np) dCTa dfA 



8 128 0.6358 0.0127 0.91% 0.45 0.61 5 1.0000 0.9809 

8 128 0.6340 0.0086 0.63% 0.46 0.71 5 0.9809 1.0000 

8 128 0.6326 0.0062 0.39% 0.40 0.77 5 0.9490 0.9910 

32 128 0.6313 0.0020 0.20% 0.62 0.54 3 0.4401 0.4357 

32 128 0.6330 0.0024 0.46% 1.20 0.77 3 0.4507 0.4630 



0.9490 0.4401 0.4507 

0.9910 0.4357 0.4630 

1.0000 0.4363 0.4639 

0.4363 1.0000 0.9267 

0.4639 0.9267 1.0000 



I'plain CTuncorr 0.6334 0.0038 0.52% 0.85 

o-corr 0.0067 0.52% 0.49 

;?orr CTuncorr 0.6322 0.0015 0.33% 1.35 

(Tcorr 0.0024 0.33% 0.84 

I>cov (Tcorr 0.6300 0.0017 -0.01% -0.05 



1.0000 1.0000 
0.0106 0.0254 
0.2485 -1.5805 



1.0000 1.0000 1.0000 
0.0503 0.5315 0.3823 
1.6625 0.7948 -0.1253 



turned out that the maxima of the various quantities un- 
der consideration were all within the reweighting range 
of this chosen simulation point for the system sizes and 
lengths of time series at hand. 

For determining the correlation-length exponent v we 
again considered the scaling of the logarithmic magneti- 
zation derivatives dln(|m|')/d/3 for i = 1, 2 and 3 and the 
derivatices of the cumulants and U^. We find scaling 
corrections to be somewhat more pronounced than for 
the two-dimensional model for the system sizes studied 
here. For the logarithmic magnetization derivatives we 
therefore performed fits of the form pT|) including the 
correction term on the full range 8 < < 128, where 
the resulting values of the effective correction exponent 
w were w = 0.57(63) {i = 1), w = 0.69(56) {i = 2) 
and w = 0.80(52) {i = 3), respectively. For the cumu- 
lants U2 and U4, on the other hand, corrections were 
too small to be fitted reliably with our data, such that 
they were effectively taken into account by dropping the 
small lattice sizes instead, while using fits of the form 
([50)1 with Ui^c = fixed. The corresponding fit data are 
collected in Table [V] The estimated standard deviations 
of the individual estimates are again found to be very 
heterogeneous, but the correlations between the different 
estimates are somewhat smaller than in two dimensions, 
in particular between the magnetization derivatives and 
the cumulants, cf. Table IVl Comparing to the case of fits 
without corrections, it is seen that this latter effect is par- 
tially due to the use of two different fit forms for the two 
types of quantities. (The fits for U2 and U4 also include a 
reduced range of lattice sizes which could lead to a decor- 
relation, but this effect is found to be much less important 
than the difference in the fit forms.) Considering the av- 
erages of individual estimates, as a result of these smaller 



correlations the underestimation of statistical errors in 
the naive approach as well as the reduction of variance 
through the optimized estimator ([22| is somewhat less 
dramatic than for the two-dimensional model, but the 
qualitative behavior appears to be very much the same. 
As our final estimate we quote ly = 0.6300(17), very 
well in agreement with the reference value v = 0.6301(4) 
taken from a survey of recent literature estimates com- 
piled in Ref. [44 1 . 



In a second step we determined the transition coupling 
from fits of the functional form (|28|) to the maxima of the 
quantities listed in Table [TTl As for the ly fits, however, 
the inclusion of an effective correction term as indicated 
in Eq. (^5)) turned out to be necessary for a faithful de- 
scription of the scaling data. The plain, error-weighted 
and covariancc-weighted averages of the corresponding 
estimates arc listed in the first two data columns of Ta- 
ble IVII together with their standard deviations, the re- 
sults being consistent with the reference value. We also 
tried non-linear three-parameter fits of the form ([^5]) to 
the data, determining v and /3c simultaneously. For this 
case, the precision of the data is not high enough to re- 
liably include corrections to scaling. Still, the improved 
results are well consistent with the reference values of 
Refs. Hlil, cf. the middle columns of Table IVll 



Finally, we also considered the scaling dimensions 
and Xt ■ For the magnetic scaling dimension, we find that 
the determinations from Xh = P /v and s/i = 3/2 — ^ /2v 
are only very weakly correlated, such that the error- 
weighted and covariance-weighted averages are very sim- 
ilar, see the right hand side of Table IVII Larger corre- 
lations are present again between the different estimates 
of the energetic scaling dimension xt from the various 
estimates of v via Xt = S — l/v and the scaling of the 
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TABLE VI: DifTerent averages for the 3D Ising model and the associated standard deviations for the transition couphng /3c 
from fits of the form (|28|l with v — 0.6301 fixed, from non-linear three-parameter fits of the form (|28[) yielding v and Pc 
simultaneously, and for the magnetic and energetic scaling dimensions according to Eq. (1341). The reference values for Xh and 
Xt have been computed from the values /3 — 0.3265(3) and u = 0.6301(4) taken from Ref. [4J| via Eq. (|34|l . 







Eq. (dlj. 


u = 0.6301 
a 


V 


a 


Eq. (HIJ 


a 




Eq. ( 

a 


131 




^plain 


•^uncorr 


0.22165681 


0.00000108 


0.6020 


0.0105 


0.2216530 


0.0000025 


0.51364 


0.00401 


1.4137 


0.0138 




corr 




0.00000170 




0.0150 




0.0000032 




0.00435 




0.0184 


^crr 


^uncorr 


0.22165741 


0.00000059 


0.6247 


0.0062 


0.2216550 


0.0000008 


0.51489 


0.00381 


1.4180 


0.0038 




^ corr 




0.00000114 




0.0077 




0.0000016 




0.00413 




0.0061 


^cov 


^ corr 


0.22165703 


0.00000085 


0.6381 


0.0044 


0.2216552 


0.0000011 


0.51516 


0.00412 


\AV1\ 


0.0043 


reference 




0.22165459 


0.00000006 


0.6301 


0.0004 


0.22165459 


0.00000006 


0.51817 


0.00058 


1.4130 


0.0010 



specific heat via Xt = 3/2 — ajlv^ leading to a consid- 
erable improvement in precision of the optimal average 
over the plain and error- weighting schemes. The results 
for both scaling dimensions are well compatible with the 
values Xh = 0.51817(58) and xt = 1.4130(10) extracted 
from the reference values of Ref. [31 . 



VI. CONCLUSIONS 

Time series data from Markov chain Monte Carlo sim- 
ulations are usually analyzed in a variety of ways to ex- 
tract estimates for the parameters of interest such as, 
e.g., critical exponents, transition temperatures, latent 
heats etc. As long as at least some of these estimates 
are based on the same simulation data, a certain degree 
of cross correlations between estimators is unavoidable. 
We have shown for the case of a finite-size scaling anal- 
ysis of the ferromagnetic nearest-neighbor Ising model 
on square and cubic lattices that more often than not, 
such correlations are very strong, with correlation coeffi- 
cients well above 0.8. While such correlations, although 
their existence is rather obvious, have been traditionally 
mostly neglected even in high-precision numerical simu- 
lation studies, it was shown here that their presence is of 
importance at different steps of the process of data anal- 
ysis, and neglecting them leads to systematically wrong 
estimates of statistical fluctuations as well as non-optimal 
combination of single estimates into final averages. 

As far as the general statistical analysis of simulation 
data is concerned, it has been discussed that traditional 
prescriptions such as error propagation have their short- 
comings, in particular as soon as non-parametric steps 
such as the determination of a maximum via reweighting 
or fitting procedures come into play. These problems are 
circumvented by resorting to the class of non-parametric 
resampling schemes, of which we have discussed the jack- 
knife technique as a conceptually and practically very 
simple representative. Using this technique, we have 
outlined a very general framework of data analysis for 
MCMC simulations consisting of (a) a transformation of 
the original set of time series into an auxiliary set of 
"binned" series, where successive samples are approxi- 



mately uncorrelated in time and (b) a general jackknif- 
ing framework, where the required steps of computing a 
parameter estimate — possibly including reweighting or 
fitting procedures etc. — are performed on the full un- 
derlying time series apart from a small window cut out 
from the data stream allowing for a reliable and robust es- 
timate of variances and covariances as well as bias effects 
without any non-stochastic approximations. While this 
technique of data analysis is not new, we feel that it still 
has not found the widespread use it deserves and hope 
that the gentle and detailed introduction given above will 
contribute to a broader adoption of this approach. 

A particular example of where the presence of cross 
correlations comes into play occurs when taking averages 
of different estimates for a parameter from the same data 
base. Neglecting correlations there leads to (a) systemat- 
ically wrong, most often too small, estimates of statisti- 
cal errors of the resulting averages and (b) a sub-optimal 
weighting of individual values in the average leading to 
larger-than-necessary variances. Correct variances can 
be estimated straightforwardly from the jackknifing ap- 
proach, while optimal weighting involves knowledge of 
the covariance matrix which is a natural byproduct of 
the jackknife technique as well. We have discussed these 
concepts in some detail for the case of a finite-size scal- 
ing analysis of the critical points of the 2D and 3D Ising 
models. It is seen there that the plain and error-weighted 
averages most oftenly used in fact can have larger fluc- 
tuations than the most precise single estimates entering 
them, but this flaw is not being detected by the con- 
ventional analysis due to the generic underestimation of 
variances. On the contrary, by using the truly optimal 
weighting of individual estimates an often substantial 
reduction of statistical fiuctuations as compared to the 
error-weighting scheme can be achieved. For some of the 
considered examples, a threefold reduction in standard 
deviation, corresponding to saving an about tenfold in- 
crease in computer time necessary to achieve the same 
result with the conventional analysis, can be achieved 
with essentially no computational overhead. In view of 
these results, heuristic rules such as, e.g., taking an error- 
weighted average using the smallest single standard devi- 
ation as an error estimate are clearly found to be inade- 
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quatc. We therefore see only two statistically acceptable 
ways of dealing with the existence of several estimates for 
the same quantity: (a) select the single most precise esti- 
mate and discard the rest or (b) combine all estimates in 
a statistically optimal way taking cross correlations into 
account. Needless to say, the latter approach is generally 
preferable in that it leads to more precise results at very 
low costs. 

We suggest to use the existence of scaling relations be- 
tween the critical exponents for the case of a continuous 
phase transition to improve the precision of estimates by 
considering the scaling dimensions as the parameters of 
primary interest. Performing the corresponding analy- 
sis taking cross correlations into account, results in a set 
of critical exponents with reduced statistical fluctuations 
that fulfill the scaling relations exactly. An application 
of this type of approach initially suggested in Ref. [l^ 
for using mean-value relations such as Callen identities 
or Schwinger-Dyson equations instead of scaling relations 
has been discussed in Ref. j20j . 

While the examples discussed were specific, it should 
be clear that the method itself is rather generic, and 
should apply to all data sets generated from MCMC sim- 
ulations. In particular, it is easy to envisage applications 
in the theory of critical phenomena, reaching from classi- 
cal statistical mechanics [i^ over soft matter physics Q 
to quantum phase transitions @, or for studying first- 
order phase transitions [l^l- The range of applications 
is not restricted to MCMC simulations, however, but 
applies with little or no modifications to other random 
sampling problems, such as, e.g. stochastic ground-state 
computations [i^, |4^ or the samplin g of polymer config- 
urations with chain-growth methods [50|, l5l| ■ 
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Appendix A: Optimal average of two correlated 
variables 



Consider a general average of two random variables xi 
and X2 [17| . 



X = KXi + (1 — k)x2, 



where < k, < 1. According to Eq. ((2T|). the variance of 
X is 

Cr^(x) = K^CTi + 2k(1 - K)pCri(T2 + (1 - k)^0-2, (A1) 

where cr^ and cr| are the variances of xi and X2; respec- 
tively, and p denotes the correlation coefficient of xi and 
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FIG. 4: (Color online) Generic form of the minimal variance 
cy'^(x) of Eq. (|A3|) as a function of the correlation coefRcient 



X2, p ~ ri2/o'icr2. Eq. (jAip is a quadratic form in k, 
which has a minimum as long as 

'^1 + ^ 2p(Tl(T2 > 0, 

which is almost always fulfilled since — 1 < p < 1: 
^1 + (^2 ~ 2p(Ti(T2 > (cTi — 0-2)^ > 0. 

Equality holds only for di = (T2 = cr and p = 1 . in which 
case any choice of k yields the same variance <j'^{x) = . 
In all other cases, the optimal weights are 



1 - K = 



l/g| - p/uiU2 

\/al + l/cr| - 2plaia2 



(A2) 



l/(jl + l/al-2p/{aia2y 
and the resulting variance of the average is 

l/al + l/al-2p/{cj^a2y 
A number of observations are immediate 



(A3) 



(i) For the uncorrelated case p = 0, one arrives back 
at the error- weighted average of Eqs. (|25p and (|26p . 

(ii) In the correlated case, and for fixed variances a\ 
and (t|, the variance cr^(i) smoothly depends on 
the correlation coefficient p. It has maxima at 
<Ji I <J2 and (72 /(Ti , only one of which is in the range 
IpI ^ 1- Notably, the relevant maximum is always 
at non-negative values of p. 

(iii) For p = ±1, the variance vanishes identically, apart 
from the singular case cti = (T2 and p = 1. 

The generic form of cr'^{x) as a function of p is depicted 
in Fig. 21 In the presence of moderate correlations, there- 
fore, anti- correlations are preferable over correlations in 
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terms of reducing the variance of the average. Note 
that the result (|A3p is different from that of Eq. (8) 
in Ref. [lOl, since the definition of correlation coefficient 
used there is different from that in our situation of taking 
an average. Instead of measuring the correlation between 
xi and X2, their definition refers to the correlation of xi 
and X2 — xi. 

The weights k and 1 — k of Eq. (jA2p are not restricted 
to be between zero and one. It is easy to see that for 
K > 1 or K < 0, the average x is in fact outside of 
the bracket [min(a;i, X2), max(xi, .T2)]. This seemingly 
paradoxical effect is easily understood from the optimal 
weights derived here. From Eq. (jA2p one reads off that 
the weights k and 1 — k leave the range 0<k,1 — k<1 
as soon as p > oxjai resp. p > a2/o'i, depending on 
whether ai < ai or a^. < cri, that is, only for strong pos- 
itive correlations to the right of the maximum in Fig. |4l 
Thus, if the smaller of a;i and X2 has the smaller vari- 
ance (and both are strongly correlated), the average is 
below both values. If the larger value has the smaller 
variance, the optimal average is above both values. The 
asymmetry comes here from the difference in variance. 



To understand this intuitively, assume for instance that 
x\ < X2 and cri < a2 with strong positive correlations 
p > iJi/(j2. It is most likely, then, that xi and X2 devi- 
ate in the same direction from the true mean (x) . Since 
f 1 < 172, the deviation of xi should be generically smaller 
than that of X2- For xi < X2, however, this is only pos- 
sible if (x) < xi < X2- This is illustrated in Fig. [S) 



X2 



Xl 

{X) 

FIG. 5: For strong positive correlations, i.e., for p > 0-1/02 in 
the case o\ < 02, the most likely location of the true expec- 
tation (a;) is outside of the bracket [min(xi, X2), max(a:i, 3:2)]. 
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